#!/usr/bin/env python3
# -*- coding: utf-8 -*-

# -------------------------------------------------------------------------------

# This file is part of code_saturne, a general-purpose CFD tool.
#
# Copyright (C) 1998-2024 EDF S.A.
#
# This program is free software; you can redistribute it and/or modify it under
# the terms of the GNU General Public License as published by the Free Software
# Foundation; either version 2 of the License, or (at your option) any later
# version.
#
# This program is distributed in the hope that it will be useful, but WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more
# details.
#
# You should have received a copy of the GNU General Public License along with
# this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
# Street, Fifth Floor, Boston, MA 02110-1301, USA.

# -------------------------------------------------------------------------------

import os
import re

from code_saturne.base.cs_math_parser import cs_math_parser

from code_saturne.model.NotebookModel import NotebookModel
from code_saturne.model.TimeTablesModel import TimeTablesModel
from code_saturne.model.SolutionDomainModel import getRunType

# ===============================================================================
# Code block templates
# ===============================================================================

# -------------------------------------------------------------------------------

_file_header = \
"""/*----------------------------------------------------------------------------*/
/*
  This file is generated by code_saturne, a general-purpose CFD tool.
*/
/*----------------------------------------------------------------------------*/

#include "base/cs_defs.h"

/*----------------------------------------------------------------------------
 * Standard C library headers
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <math.h>

#if defined(HAVE_MPI)
#include <mpi.h>
#endif

/*----------------------------------------------------------------------------
 *  Local headers
 *----------------------------------------------------------------------------*/

#include "cs_headers.h"
"""

_file_header2 = \
"""#include "nc_phase.h"

"""

_file_header3 = \
"""
/*----------------------------------------------------------------------------*/

BEGIN_C_DECLS

/*----------------------------------------------------------------------------*/

"""

_private_func_header = \
"""
/* =============================================================================
 * Private functions
 * ===========================================================================*/

"""

_public_func_header = \
"""
/* =============================================================================
 * Public functions
 * ===========================================================================*/

"""

_function_close = \
"""}

"""

_file_footer = \
"""}

/*----------------------------------------------------------------------------*/

END_C_DECLS

"""

_function_getter_header = { \
'ibm':"""cs_cutcell_func_t *
cs_meg_ibm_func_by_name(const char *object_name)
{
""",
'ibm_vol':"""cs_ibm_volume_func_t *
cs_meg_ibm_volume_func_by_name(const char *object_name,
                               const char *var_name)
{
""",
'ibm_fsi':"""cs_ibm_fsi_func_t *
cs_meg_ibm_fsi_func_by_name(const char *object_name,
                            const char *var_name)
{
"""
}

_internal_function_name = { \
'ibm':'_meg_immersed_boundaries',
'ibm_vol':'_meg_immersed_boundaries_volume',
'ibm_fsi':'_meg_immersed_boundaries_fsi'
}

_internal_function_prefix = { \
'ibm':"""/*----------------------------------------------------------------------------*/

static int
{_name}
(
  [[maybe_unused]] const cs_lnum_t c_id,
  [[maybe_unused]] const cs_real_t xyz[3],
  [[maybe_unused]] const cs_real_t time,
  [[maybe_unused]] const int       num_object
)
{_sep_open}
""",
'ibm_vol':"""/*----------------------------------------------------------------------------*/

static void
{_name}
(
  [[maybe_unused]] const cs_lnum_t c_id,
  [[maybe_unused]] const cs_real_t xyz[3],
  [[maybe_unused]] const cs_real_t time,
                   cs_real_t      *retvals
)
{_sep_open}
""",
'ibm_fsi':"""/*----------------------------------------------------------------------------*/

static void
{_name}
(
  cs_real_t      *retvals
)
{_sep_open}
"""
}

_function_header = { \
'vol':"""void
cs_meg_volume_function(const char      *zone_name,
                       [[maybe_unused]] const cs_lnum_t  n_elts,
                       [[maybe_unused]] const cs_lnum_t *elt_ids,
                       [[maybe_unused]] const cs_real_t  xyz[][3],
                       const char      *fields_names,
                       cs_real_t       *fvals[])
{
""",
'bnd':"""void
cs_meg_boundary_function(const char       *zone_name,
                         [[maybe_unused]] const cs_lnum_t   n_elts,
                         [[maybe_unused]] const cs_lnum_t  *elt_ids,
                         [[maybe_unused]] const cs_real_t   xyz[][3],
                         const char       *field_name,
                         [[maybe_unused]] const char       *condition,
                         cs_real_t        *retvals)
{
""",
'src':"""void
cs_meg_source_terms(const char       *zone_name,
                    [[maybe_unused]] const cs_lnum_t   n_elts,
                    [[maybe_unused]] const cs_lnum_t  *elt_ids,
                    [[maybe_unused]] const cs_real_t   xyz[][3],
                    const char       *name,
                    [[maybe_unused]] const char       *source_type,
                    cs_real_t        *retvals)
{
""",
'ini':"""void
cs_meg_initialization(const char      *zone_name,
                      [[maybe_unused]] const cs_lnum_t  n_elts,
                      [[maybe_unused]] const cs_lnum_t *elt_ids,
                      [[maybe_unused]] const cs_real_t  xyz[][3],
                      const char      *field_name,
                      cs_real_t       *retvals)
{
""",
'fsi':"""void
cs_meg_fsi_struct(const char       *object_type,
                  const char       *name,
                  const cs_real_t   fluid_f[],
                  cs_real_t         val[])
{
""",
'pfl':"""void
cs_meg_post_profiles(const char       *name,
                     int               n_coords,
                     cs_real_t         coords[][3])
{
""",
'pwa':"""void
cs_meg_post_activate(void)
{
""",
'pca':"""void
cs_meg_post_calculator(const char       *name,
                       [[maybe_unused]] const cs_lnum_t   n_elts,
                       [[maybe_unused]] const cs_lnum_t  *elt_ids,
                       [[maybe_unused]] const cs_real_t   xyz[][3],
                       cs_real_t        *retvals)
{
"""
}

_function_names = {'vol': 'cs_meg_volume_function.cxx',
                   'bnd': 'cs_meg_boundary_function.cxx',
                   'src': 'cs_meg_source_terms.cxx',
                   'ini': 'cs_meg_initialization.cxx',
                   'ibm': 'cs_meg_immersed_boundaries.cxx',
                   'ibm_vol': 'cs_meg_immersed_boundaries_volume.cxx',
                   'ibm_fsi': 'cs_meg_immersed_boundaries_fsi.cxx',
                   'fsi': 'cs_meg_fsi_struct.cxx',
                   'pfl': 'cs_meg_post_profile.cxx',
                   'pwa': 'cs_meg_post_output.cxx',
                   'pca': 'cs_meg_post_calculator.cxx'}

_block_comments = {'vol': 'User defined formula for variable(s) %s over zone %s',
                   'bnd': 'User defined formula for "%s" over BC=%s',
                   'src': 'User defined source term for %s over zone %s',
                   'ini': 'User defined initialization for variable %s over zone %s',
                   'ibm': 'User defined explicit formula of %s indicator for object %s',
                   'ibm_vol': 'User defined physical property for variable %s for IBM object %s',
                   'ibm_fsi': 'User defined FSI coupling structure %s for IBM object %s',
                   'fsi': 'User defined FSI coupling structure %s for zone %s',
                   'pfl': 'User-defined %s for profile %s',
                   'pwa': 'User defined %s for writer %s',
                   'pca': 'User defined postprocessing formula for "%s"'}

_func_short_to_long = {'vol': 'volume zone',
                       'bnd': 'boundary',
                       'src': 'source term',
                       'ini': 'initialization',
                       'ibm': 'Immersed boundaries explicit function',
                       'ibm_vol': 'Immersed volume zone',
                       'ibm_fsi': 'Immersed FSI',
                       'fsi': 'Mechanicaly-coupled structures',
                       'pfl': 'Profile coordinates',
                       'pwa': 'Writer activation',
                       'pca': 'Postprocessing calculator'}

# -------------------------------------------------------------------------------

_pkg_fluid_prop_dict = {}
_pkg_fluid_prop_dict['code_saturne'] = {'rho0':'ro0',
                                        'mu0':'viscl0',
                                        'p0':'p0',
                                        't0':'t0',
                                        'cp0':'cp0',
                                        'lambda0':'lambda0',
                                        'viscv0':'viscv0'}

_pkg_fluid_prop_dict['neptune_cfd'] = {'rho0':'ro0',
                                       'mu0':'viscl0',
                                       'cp0':'cp0',
                                       'lambda0':'lambda0'}

_pkg_glob_struct = {'code_saturne':'cs_glob_fluid_properties',
                    'neptune_cfd':'nc_phases->p_ini[PHASE_ID]'}

# ---------------------------------------------------------------------------

_base_tokens = {'dt':'const cs_real_t dt = cs_glob_time_step->dt[0];',
                't':'const cs_real_t t = cs_glob_time_step->t_cur;',
                'iter':'const int iter = cs_glob_time_step->nt_cur;',
                'volume':'const cs_real_t volume = cs_volume_zone_by_name(zone_name)->measure;',
                'fluid_volume':'const cs_real_t fluid_volume = cs_volume_zone_by_name(zone_name)->f_measure;',
                'surface':'const cs_real_t surface = cs_boundary_zone_by_name(zone_name)->measure;',
                'fluid_surface':'const cs_real_t fluid_surface = cs_boundary_zone_by_name(zone_name)->f_measure;',
                'pi':'const cs_real_t pi = cs_math_pi;',
                'uref':'const cs_real_t uref = cs_glob_turb_ref_values->uref;',
                'almax':'const cs_real_t almax = cs_glob_turb_ref_values->almax;',
                'gx':'const cs_real_t gx = cs_glob_physical_constants->gravity[0];',
                'gy':'const cs_real_t gy = cs_glob_physical_constants->gravity[1];',
                'gz':'const cs_real_t gz = cs_glob_physical_constants->gravity[2];'}

# ---------------------------------------------------------------------------

def _error_and_exit(msg):
    """
    Print error message to stderr and exit.
    """

    import sys

    print(msg, file=sys.stderr)
    sys.exit(1)

# ---------------------------------------------------------------------------

def build_internal_function_header(func_type, suffix):
    """
    Build a header for an internal function
    """

    _int_f_name = '_'.join([_internal_function_name[func_type], suffix])

    return _int_f_name


# ---------------------------------------------------------------------------

def parse_gui_expression(expression,
                         req,
                         known_symbols,
                         func_type,
                         glob_tokens,
                         loop_tokens,
                         need_for_loop = False,
                         indent_decl = 2,
                         indent_main = 3):
    """
    Parse GUI expression.
    """

    usr_code = ''
    usr_defs = ''
    if_loop = False

    tab = '  '
    ntabs = indent_main
    if not need_for_loop:
        ntabs -= 1


    parser = cs_math_parser()

    expr_list, expr_user = parser.parse_expression(expression,
                                                   req,
                                                   known_symbols,
                                                   func_type,
                                                   glob_tokens,
                                                   loop_tokens,
                                                   need_for_loop)

    for exp in expr_user:
        usr_defs += indent_decl*tab + exp

    for exp in expr_list:
        usr_code += ntabs*tab + exp

    return usr_code, usr_defs

# ===============================================================================
# Utility functions
# ===============================================================================

def break_expression(exp):

    expression_lines = []

    for line in exp.split('\n'):
        line_comp = []
        for elt in re.split('=|\+|-|\*|\/|\(|\)|;|,|\^|<|>|\&\&|\|\|',line):
            if elt != '':
                line_comp.append(elt.strip())

        expression_lines.append(line_comp)

    return expression_lines

# ===============================================================================
# Main class
# ===============================================================================

class meg_to_c_interpreter:

    # ---------------------------------------------------------------------------

    def __init__(self,
                 case,
                 create_functions=True,
                 module_name=None,
                 wdir=None):

        self.case = case
        self.wdir = wdir

        if module_name:
            self.module_name = module_name
        else:
            self.module_name = case.module_name()

        if not self.wdir:
            data_path = os.path.join(case['case_path'], 'DATA')
        else:
            data_path = self.wdir

        self.tmp_path = os.path.join(data_path, 'tmp')

        # function name to file name dictionary
        self.funcs = {'vol': {},
                      'bnd': {},
                      'src': {},
                      'ini': {},
                      'ibm': {},
                      'ibm_vol': {},
                      'ibm_fsi': {},
                      'fsi': {},
                      'pfl': {},
                      'pwa': {},
                      'pca': {}}

        self.code_to_write = ""

        nb = NotebookModel(self.case)
        self.notebook = {}
        for (nme, val) in nb.getNotebookList():
            self.notebook[nme] = str(val)

        TTM = TimeTablesModel(self.case)
        self.time_tables = {}
        for i_t, tab_name in enumerate(TTM.getTableNamesList()):
            for _h in TTM.getTableHeadersList(i_t):
                _k = "{}[{}]".format(tab_name, _h)
                self.time_tables[_k] = ['{}_{}'.format(tab_name,_h),
                                        'CS_TIME_TABLE("{}","{}")'.format(tab_name, _h)]

        if create_functions and getRunType(self.case) == 'standard':

            from code_saturne.model.LocalizationModel import LocalizationModel

            vlm = LocalizationModel('VolumicZone', self.case)
            vol_zones = vlm.getZones()

            # Volume code
            self.generate_volume_code(vol_zones)

            # Boundary code
            self.generate_boundary_code()

            # Source terms code
            self.generate_source_terms_code(vol_zones)

            # Initialization function
            self.generate_initialize_code(vol_zones)

            # Immersed boundaries function
            self.generate_immersed_boundaries_code()

            # ALE/FSI internal coupling function
            self.generate_fsi_ic_code()

            # Writer profiles code
            self.generate_post_profile_code()

            # Writer activation
            self.generate_writer_activation_code()

            # Postprocessing calculator
            self.generate_post_calculator_code()

    # ---------------------------------------------------------------------------

    def update_block_expression(self, func_type, key, new_exp):

        self.funcs[func_type][key]['exp']   = new_exp
        self.funcs[func_type][key]['lines'] = break_expression(new_exp)

    # ---------------------------------------------------------------------------

    def init_block(self,
                   ftype,
                   zone_name,
                   name,
                   expression,
                   required,
                   symbols,
                   known_fields,
                   condition = None,
                   source_type = None,
                   element_type= "center"):

        # Creating a unique function name based on the zone and variable name
        fkey = '::'.join([zone_name, name])
        if fkey in self.funcs[ftype].keys():
            msg = 'Formula for "%s" in %s %s was already defined:\n %s' \
                    % (name, _func_short_to_long[ftype], zone_name,
                       self.funcs[ftype][fkey]['exp'])
            raise Exception(msg)


        self.funcs[ftype][fkey] = {'exp': expression,
                                   'req': required,
                                   'sym': symbols,
                                   'knf': known_fields,
                                   'cnd': condition,
                                   'tpe': source_type,
                                   'elt': element_type}
        if self.funcs[ftype][fkey]['elt'] not in ('center', 'vertex'):
            self.funcs[ftype][fkey]['elt'] = 'center'

        self.funcs[ftype][fkey]['lines'] = break_expression(expression)

    # ---------------------------------------------------------------------------

    def write_cell_block(self, func_key):

        func_params = self.funcs['vol'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        zone, name = func_key.split('::')
        exp_lines_comp = func_params['lines']

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']
        need_coords = False

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[c_id][%s];' % (kc, str(ic))

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # fluid properties
        for kp in _pkg_fluid_prop_dict[self.module_name].keys():
            if len(name.split("_")) > 1:
                try:
                    phase_id = int(name.split('_')[-1])-1
                except:
                    phase_id = -1
            else:
                phase_id = -1
            gs = _pkg_glob_struct[self.module_name].replace('PHASE_ID',
                                                            str(phase_id))
            pn = _pkg_fluid_prop_dict[self.module_name][kp]
            glob_tokens[kp] = 'const cs_real_t %s = %s->%s;' %(kp, gs, pn)

        if name[-12:] == '_diffusivity':
            name_ref = name + '_ref'
            for s in symbols:
                if s[0] == name_ref:
                    base_name = name[:-12]
                    glob_tokens[name_ref] = 'const cs_real_t %s\n' % (name_ref)
                    glob_tokens[name_ref] += '      = cs_field_get_key_double(cs_field_by_name("%s"), cs_field_key_id("diffusivity_ref"));' % (base_name)

        # Fields
        label_not_name = ['Additional scalar', 'Thermal scalar', 'Pressure']
        for f in known_fields:
            # Get label, name and dimension
            try:
                (fl, fn, fdim) = f
                if fdim < 1:
                    fdim = 1
            except:
                (fl, fn) = f
                fdim = 1

            for lnn in label_not_name:
                if lnn in fn:
                    fn = fl

            glob_tokens[fl] = \
                 'const cs_real_t *%s_vals = cs_field_by_name("%s")->val;' \
                 % (fl, fn)

            if fdim == 1:
                loop_tokens[fl] = 'const cs_real_t %s = %s_vals[c_id];' \
                     % (fl, fl)
            else:
                loop_tokens[fl] = 'const cs_real_t *%s = %s_vals + %d*c_id;' \
                     % (fl, fl, fdim)

        # ------------------------

        for s in required:
            known_symbols.append(s);

        known_symbols.append('#')

        ntabs += 1

        # Parse the user expression
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'vol',
                                          glob_tokens,
                                          loop_tokens,
                                          need_for_loop=True)
        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        nsplit = name.split('+')
        usr_blck = tab + 'if (strcmp(fields_names, "%s") == 0 &&\n' % (name)
        usr_blck += tab + '    strcmp(zone_name, "%s") == 0) {\n' % (zone)

        usr_blck += usr_defs

        usr_blck += 2*tab + 'for (cs_lnum_t e_id = 0; e_id < n_elts; e_id++) {\n'
        usr_blck += 3*tab + 'cs_lnum_t c_id = (elt_ids != nullptr) ? elt_ids[e_id] : e_id;\n'

        usr_blck += usr_code

        usr_blck += 2*tab + '}\n'
        usr_blck += tab + '}\n'

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_bnd_block(self, func_key):

        if func_key not in self.funcs['bnd'].keys():
            return None

        func_params = self.funcs['bnd'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']
        cname        = func_params['cnd']
        element_type = func_params['elt']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        exp_lines_comp = func_params['lines']

        zone, field_name = func_key.split('::')

        # Check if for loop is needed
        need_for_loop = True
        if cname == 'flow1_formula' or cname == 'flow2_formula':
            need_for_loop = False

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        for req in required:
            known_symbols.append(req)

        coords = ['x', 'y', 'z']
        need_coords = False

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[b_e_id][%s];' \
                    % (kc, str(ic))

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Time table variables
        for ktt in self.time_tables.keys():
            _kvals = self.time_tables[ktt]
            glob_tokens[ktt] = \
                'const cs_real_t {var} = {func};'.format(var=_kvals[0], \
                                                          func=_kvals[1])

        # Fields
        for f in known_fields:
            glob_tokens[f[0]] = \
            'const cs_real_t *%s_vals = cs_field_by_name("%s")->val;' % (f[0], f[1])
            loop_tokens[f[0]] = \
            'const cs_real_t %s = %s_vals[b_e_id];' % (f[0], f[0])
        # ------------------------

        if need_for_loop:
            ntabs += 1

        # Parse the user expression
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'bnd',
                                          glob_tokens,
                                          loop_tokens,
                                          need_for_loop)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        block_cond  = tab + 'if (strcmp(field_name, "%s") == 0 &&\n' % (field_name)
        block_cond += tab + '    strcmp(condition, "%s") == 0 &&\n' % (cname)
        block_cond += tab + '    strcmp(zone_name, "%s") == 0) {\n' % (zone)
        usr_blck = block_cond + '\n'

        usr_blck += usr_defs

        if need_for_loop:
            usr_blck += 2*tab + 'for (cs_lnum_t e_id = 0; e_id < n_elts; e_id++) {\n'
            usr_blck += 3*tab + 'cs_lnum_t b_e_id = (elt_ids != nullptr) ? elt_ids[e_id] : e_id;\n'

        usr_blck += usr_code

        if need_for_loop:
            usr_blck += 2*tab + '}\n'

        usr_blck += tab + '}\n'

        # Replace time table calls
        for _tt in self.time_tables.keys():
            usr_blck = usr_blck.replace(_tt, self.time_tables[_tt][0])

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_src_block(self, func_key):

        # Check if function exists:
        if func_key not in self.funcs['src'].keys():
            return

        func_params = self.funcs['src'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']

        known_fields = dict( zip([k[0] for k in func_params['knf']],
                                 [k[1] for k in func_params['knf']]))
        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        source_type  = func_params['tpe']

        zone, name = func_key.split('::')
        exp_lines_comp = func_params['lines']

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = '  '
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[c_id][%s];' % (kc, str(ic))

        # For momentum also define u,v and w:
        if source_type == "momentum_source_term":
            glob_tokens['velocity'] = \
            'const cs_real_3_t *vel = (cs_real_3_t *)CS_F_(vel)->val;'
            for i, key in enumerate(['u', 'v', 'w']):
                loop_tokens[key] = \
                'const cs_real_t %s = vel[c_id][%d];' % (key, i)

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Time table variables
        for ktt in self.time_tables.keys():
            _kvals = self.time_tables[ktt]
            glob_tokens[ktt] = \
                'const cs_real_t {var} = {func};'.format(var=_kvals[0], \
                                                         func=_kvals[1])

        for f in known_fields.keys():
            knf_name = known_fields[f]
            glob_tokens[f] = \
            'const cs_real_t *%s_vals = cs_field_by_name("%s")->val;' \
            % (f, knf_name)

            loop_tokens[f] = 'const cs_real_t %s = %s_vals[c_id];' % (f, f)

        # ------------------------

        for r in required:
            known_symbols.append(r)

        # Parse the user expression
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'src',
                                          glob_tokens,
                                          loop_tokens,
                                          need_for_loop=True)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        block_cond  = tab + 'if (strcmp(zone_name, "%s") == 0 &&\n' % (zone)
        block_cond += tab + '    strcmp(name, "%s") == 0 && \n' % (name)
        block_cond += tab + '    strcmp(source_type, "%s") == 0) {\n' % (source_type)
        usr_blck = block_cond + '\n'

        usr_blck += usr_defs

        usr_blck += 2*tab + 'for (cs_lnum_t e_id = 0; e_id < n_elts; e_id++) {\n'
        usr_blck += 3*tab + '[[maybe_unused]] cs_lnum_t c_id = (elt_ids != nullptr) ? elt_ids[e_id] : e_id;\n'

        usr_blck += usr_code

        usr_blck += 2*tab + '}\n'
        usr_blck += tab + '}\n'

        # Replace time table calls
        for _tt in self.time_tables.keys():
            usr_blck = usr_blck.replace(_tt, self.time_tables[_tt][0])

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_ini_block(self, func_key):

        # Check if function exists:
        if func_key not in self.funcs['ini'].keys():
            return

        func_params = self.funcs['ini'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']

        if (len(func_params['req']) > 0):
            if type(func_params['req'][0]) == tuple:
                required = [r[0] for r in func_params['req']]
            else:
                required = func_params['req']
        else:
            required = []

        zone, name = func_key.split('::')
        exp_lines_comp = func_params['lines']

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = '  '
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[c_id][%s];' % (kc, str(ic))

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # fluid properties
        for kp in _pkg_fluid_prop_dict[self.module_name].keys():
            if len(name.split("_")) > 1:
                try:
                    phase_id = int(name.split('_')[-1])-1
                except:
                    phase_id = -1
            else:
                phase_id = -1
            gs = _pkg_glob_struct[self.module_name].replace('PHASE_ID',
                                                         str(phase_id))
            pn = _pkg_fluid_prop_dict[self.module_name][kp]
            glob_tokens[kp] = 'const cs_real_t %s = %s->%s;' %(kp, gs, pn)

        # known fields
        for f in known_fields:
            glob_tokens[f[0]] = \
            'const cs_real_t *%s_vals = cs_field_by_name("%s")->val;' % (f[0], f[1])
            loop_tokens[f[0]] = \
            'const cs_real_t %s = %s_vals[c_id];' % (f[0], f[0])

        # ------------------------

        for r in required:
            known_symbols.append(r)

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'ini',
                                          glob_tokens,
                                          loop_tokens,
                                          need_for_loop=True)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        block_cond  = tab + 'if (strcmp(zone_name, "%s") == 0 &&\n' % (zone)
        block_cond += tab + '    strcmp(field_name, "%s") == 0) {\n' % (name)
        usr_blck = block_cond + '\n'

        usr_blck += usr_defs

        usr_blck += 2*tab + 'for (cs_lnum_t e_id = 0; e_id < n_elts; e_id++) {\n'
        usr_blck += 3*tab + '[[maybe_unused]] cs_lnum_t c_id = (elt_ids != nullptr) ? elt_ids[e_id] : e_id;\n'

        usr_blck += usr_code

        usr_blck += 2*tab + '}\n'
        usr_blck += tab + '}\n'

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_ibm_block(self, func_key):

        func_params = self.funcs['ibm'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        object_name, name = func_key.split('::')
        exp_lines_comp = func_params['lines']

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[%s];' % (kc, str(ic))

        # For FSI objects:
        #Center of gravity
        cog_fsi = ['cog_fsi_x', 'cog_fsi_y', 'cog_fsi_z']

        for ic, kc in enumerate(cog_fsi):
            loop_tokens[kc] = 'const cs_real_t {_kc} = '.format(_kc=kc) + \
                'nc_fsi_object->cdg[num_object][{_ic}];'.format(_ic=ic)

        #rotation matrix (rot_m)
        loop_tokens['rot_fsi'] = 'cs_real_3_t *rot_fsi = '\
            'nc_fsi_object->rot_m[num_object];'

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # ------------------------

        for s in required:
            known_symbols.append(s);

        known_symbols.append('#')

        if_loop = False

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'ibm',
                                          glob_tokens,
                                          loop_tokens,
                                          indent_decl=1,
                                          indent_main=2)

        usr_blck = tab + 'int retval = 0;\n'

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]
        if usr_defs != '':
            usr_blck += usr_defs + '\n'
        usr_blck += usr_code

        usr_blck += '\n'
        usr_blck += tab + 'return retval;\n'

        return usr_blck

    #---------------------------------------------------------------------------

    def write_ibm_vol_block(self, func_key):

        func_params = self.funcs['ibm_vol'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        object_name, name = func_key.split('::')
        exp_lines_comp = func_params['lines']

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[%s];' % (kc, str(ic))

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Temperature for rho, cp, lambda
        glob_tokens["temperature"] = 'const cs_real_t temperature = nc_poro_var->temp_sol[c_id];'

        # ------------------------

        for s in required:
            known_symbols.append(s);

        known_symbols.append('#')

        ntabs += 1
        #if_loop = False

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'ibm_vol',
                                          glob_tokens,
                                          loop_tokens,
                                          indent_decl=1,
                                          indent_main=2)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]
        if usr_defs != '':
            usr_blck += usr_defs + '\n'

        usr_blck += usr_code

        return usr_blck


    #---------------------------------------------------------------------------

    def write_ibm_fsi_block(self, func_key):

        func_params = self.funcs['ibm_fsi'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        object_name, mat_name = func_key.split('::')
        exp_lines_comp = func_params['lines']

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']

        # ------------------------

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t %s = xyz[%s];' % (kc, str(ic))

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # ------------------------

        for s in required:
            known_symbols.append(s);

        known_symbols.append('#')

        ntabs += 1
        #if_loop = False

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'ibm_fsi',
                                          glob_tokens,
                                          loop_tokens,
                                          need_for_loop=True)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]
        if usr_defs != '':
            usr_blck += usr_defs + '\n'

        set_blck = ''

        if mat_name == "porous_inertia":
            usr_blck += tab*2 + 'cs_real_t i11 = 0., i12 = 0., i13 = 0.;\n'
            usr_blck += tab*2 + 'cs_real_t i22 = 0., i23 = 0., i33 = 0.;\n\n'
            set_blck += tab*2 + 'retvals[0] = i11; retvals[1] = i12; retvals[2] = i13;\n'
            set_blck += tab*2 + 'retvals[3] = i12; retvals[4] = i22; retvals[5] = i23;\n'
            set_blck += tab*2 + 'retvals[6] = i13; retvals[7] = i23; retvals[8] = i33;\n'

        usr_blck += usr_code
        usr_blck += set_blck

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_fsi_block(self, func_key):

        if func_key not in self.funcs['fsi'].keys():
            return

        func_params = self.funcs['fsi'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']
        cname        = func_params['cnd']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        exp_lines_comp = func_params['lines']

        zone, mat_name = func_key.split('::')

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        for req in required:
            known_symbols.append(req)

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'fsi',
                                          glob_tokens,
                                          loop_tokens,
                                          indent_decl=3)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        block_cond  = tab + 'if (   strcmp(object_type, "%s") == 0\n' % (mat_name)
        block_cond += tab + '    && strcmp(name, "%s") == 0) {\n' % (zone)
        usr_blck = block_cond + '\n'
        set_blck = '    }\n\n'

        if mat_name == "mass_matrix":
            usr_blck += '    cs_real_t m11 = 0, m12 = 0, m13 = 0;\n'
            usr_blck += '    cs_real_t m21 = 0, m22 = 0, m23 = 0;\n'
            usr_blck += '    cs_real_t m31 = 0, m32 = 0, m33 = 0;\n\n'
            set_blck += '    val[0] = m11; val[1] = m12; val[2] = m13;\n'
            set_blck += '    val[3] = m21; val[4] = m22; val[5] = m23;\n'
            set_blck += '    val[6] = m31; val[7] = m32; val[8] = m33;\n'
        elif mat_name == "stiffness_matrix":
            usr_blck += '    cs_real_t k11 = 0, k12 = 0, k13 = 0;\n'
            usr_blck += '    cs_real_t k21 = 0, k22 = 0, k23 = 0;\n'
            usr_blck += '    cs_real_t k31 = 0, k32 = 0, k33 = 0;\n\n'
            set_blck += '    val[0] = k11; val[1] = k12; val[2] = k13;\n'
            set_blck += '    val[3] = k21; val[4] = k22; val[5] = k23;\n'
            set_blck += '    val[6] = k31; val[7] = k32; val[8] = k33;\n'
        elif mat_name == "damping_matrix":
            usr_blck += '    cs_real_t c11 = 0, c12 = 0, c13 = 0;\n'
            usr_blck += '    cs_real_t c21 = 0, c22 = 0, c23 = 0;\n'
            usr_blck += '    cs_real_t c31 = 0, c32 = 0, c33 = 0;\n\n'
            set_blck += '    val[0] = c11; val[1] = c12; val[2] = c13;\n'
            set_blck += '    val[3] = c21; val[4] = c22; val[5] = c23;\n'
            set_blck += '    val[6] = c31; val[7] = c32; val[8] = c33;\n'
        elif mat_name == "fluid_force":
            usr_blck += '    cs_real_t fluid_fx = fluid_f[0];\n'
            usr_blck += '    cs_real_t fluid_fy = fluid_f[1];\n'
            usr_blck += '    cs_real_t fluid_fz = fluid_f[2];\n'
            usr_blck += '    cs_real_t fx = 0, fy = 0, fz = 0;\n'
            set_blck += '    val[0] = fx; val[1] = fy; val[2] = fz;\n'

        usr_blck += '    {\n'

        usr_blck += usr_defs

        usr_blck += usr_code

        usr_blck += set_blck

        usr_blck += tab + '}\n'

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_profile_coo_block(self, func_key):

        if func_key not in self.funcs['pfl'].keys():
            return

        func_params = self.funcs['pfl'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']
        cname        = func_params['cnd']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        exp_lines_comp = func_params['lines']

        name, s = func_key.split('::')

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        for req in required:
            known_symbols.append(req)

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'pwa',
                                          glob_tokens,
                                          loop_tokens,
                                          indent_decl=2)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        usr_blck =  tab + 'if (strcmp(name, "%s") == 0) {\n' % (name)
        usr_blck += usr_defs
        usr_blck += tab + '  cs_real_t x, y, z;\n\n'
        usr_blck += tab + '  for (int p_id = 0; p_id < n_coords; p_id++) {\n'
        usr_blck += tab + '    cs_real_t s = (cs_real_t)p_id / (cs_real_t)(n_coords-1);\n\n'

        usr_blck += usr_code

        usr_blck += '\n'
        usr_blck += tab + '    coords[p_id][0] = x;\n'
        usr_blck += tab + '    coords[p_id][1] = y;\n'
        usr_blck += tab + '    coords[p_id][2] = z;\n'
        usr_blck += tab*2 + '}\n'
        usr_blck += tab + '}\n'

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_writer_activation_block(self, func_key):

        if func_key not in self.funcs['pwa'].keys():
            return

        func_params = self.funcs['pwa'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']
        cname        = func_params['cnd']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        exp_lines_comp = func_params['lines']

        w_id, s = func_key.split('::')

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        for req in required:
            known_symbols.append(req)

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Parse the user expresion
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'pwa',
                                          glob_tokens,
                                          loop_tokens,
                                          indent_decl=3,
                                          indent_main=4)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        usr_blck  = tab + '{\n'
        usr_blck += tab + '  bool is_active = false;\n\n'
        usr_blck += tab*2 + '{\n'

        usr_blck += usr_defs

        usr_blck += usr_code

        usr_blck += tab*2 + '}\n'
        usr_blck += '\n'
        usr_blck += tab + '  cs_post_activate_writer_if_enabled('+w_id+', is_active);\n'
        usr_blck += tab + '}\n'

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_post_calculator_block(self, func_key):

        if func_key not in self.funcs['pca'].keys():
            return

        func_params = self.funcs['pca'][func_key]

        expression   = func_params['exp']
        symbols      = func_params['sym']
        known_fields = func_params['knf']

        if type(func_params['req'][0]) == tuple:
            required = [r[0] for r in func_params['req']]
        else:
            required = func_params['req']

        exp_lines_comp = func_params['lines']

        name = func_key.split('::')[1]

        # Get user definitions and code
        usr_defs = ''
        usr_code = ''
        usr_blck = ''

        tab   = "  "
        ntabs = 2

        known_symbols = []
        coords = ['x', 'y', 'z']
        need_coords = False

        # Deal with tokens which require a definition
        glob_tokens = {}
        loop_tokens = {}
        glob_tokens.update(_base_tokens)

        # Coordinates:
        for kc in coords:
            ic = coords.index(kc)
            loop_tokens[kc] = 'const cs_real_t {} = xyz[c_id][{}];'.format(kc, str(ic))

        # Notebook variables
        for kn in self.notebook.keys():
            glob_tokens[kn] = \
            'const cs_real_t %s = cs_notebook_parameter_value_by_name("%s");' % (kn, kn)

        # Fields
        label_not_name = ['Additional scalar', 'Thermal scalar', 'Pressure']
        for f in known_fields:
            # Get label, name and dimension
            try:
                (fl, fn, fdim) = f
                fdim = int(fdim)
                if fdim < 1:
                    fdim = 1
            except:
                (fl, fn) = f
                fdim = 1

            for lnn in label_not_name:
                if lnn in fn:
                    fn = fl

            glob_tokens[fl] = \
                 'const cs_real_t *%s_vals = cs_field_by_name("%s")->val;' \
                 % (fl, fn)

            if fdim == 1:
                loop_tokens[fl] = 'const cs_real_t %s = %s_vals[c_id];' \
                     % (fl, fl)
            else:
                loop_tokens[fl] = 'const cs_real_t *%s = %s_vals + %d*c_id;' \
                     % (fl, fl, fdim)

        # ------------------------

        for s in required:
            known_symbols.append(s);

        known_symbols.append('#')

        ntabs += 1

        # Parse the user expression
        parsed_exp = parse_gui_expression(expression,
                                          required,
                                          known_symbols,
                                          'pca',
                                          glob_tokens,
                                          loop_tokens,
                                          need_for_loop=True)

        usr_code += parsed_exp[0]
        if parsed_exp[1] != '':
            usr_defs += parsed_exp[1]

        # Write the block
        nsplit = name.split('+')
        usr_blck = tab + 'if (strcmp(name, "%s") == 0) {\n' % (name)

        usr_blck += usr_defs

        usr_blck += 2*tab + 'for (cs_lnum_t e_id = 0; e_id < n_elts; e_id++) {\n'
        usr_blck += 3*tab + '[[maybe_unused]] cs_lnum_t c_id = (elt_ids == NULL) ? e_id : elt_ids[e_id];\n'

        usr_blck += usr_code

        usr_blck += 2*tab + '}\n'
        usr_blck += tab + '}\n'

        return usr_blck

    # ---------------------------------------------------------------------------

    def write_block(self, func_type, key):

        # Check if function exists
        if key not in self.funcs[func_type].keys():
            return

        if func_type == 'vol':
            return self.write_cell_block(key)
        elif func_type == 'bnd':
            return self.write_bnd_block(key)
        elif func_type == 'src':
            return self.write_src_block(key)
        elif func_type == 'ini':
            return self.write_ini_block(key)
        elif func_type == 'ibm':
            return self.write_ibm_block(key)
        elif func_type == 'ibm_vol':
            return self.write_ibm_vol_block(key)
        elif func_type == 'ibm_fsi':
            return self.write_ibm_fsi_block(key)
        elif func_type == 'fsi':
            return self.write_fsi_block(key)
        elif func_type == 'pfl':
            return self.write_profile_coo_block(key)
        elif func_type == 'pwa':
            return self.write_writer_activation_block(key)
        elif func_type == 'pca':
            return self.write_post_calculator_block(key)
        else:
            return None

    # ---------------------------------------------------------------------------

    def generate_volume_code(self, vol_zones=None):
        # Ground water model enabled ?
        gwm = False

        from code_saturne.model.LocalizationModel import LocalizationModel
        from code_saturne.model.GroundwaterLawModel import GroundwaterLawModel

        if vol_zones == None:
            vlm = LocalizationModel('VolumicZone', self.case)
            vol_zones = vlm.getZones()

        if self.module_name == 'code_saturne':
            from code_saturne.model.FluidCharacteristicsModel \
                import FluidCharacteristicsModel

            fcm = FluidCharacteristicsModel(self.case)
            for (fk, sym) in fcm.lst:
                if fcm.getPropertyMode(fk) == 'user_law':
                    exp, req, sca, sym = fcm.getFormulaComponents(fk)
                    self.init_block('vol', 'all_cells', fk,
                                    exp, req, sym, sca)
                for zone in vol_zones:
                    zname = zone.getLabel()
                    z_id  = str(zone.getCodeNumber())
                    if zname != "all_cells" and \
                            zone.isNatureActivated('physical_properties'):
                        if fcm.getFormulaTry(fk, zone=z_id):
                            exp, req, sca, sym = \
                                fcm.getFormulaComponents(fk, zone=z_id)
                            self.init_block('vol', zname, fk,
                                            exp, req, sym, sca)

            slist = fcm.m_sca.getUserScalarNameList()
            for s in fcm.m_sca.getScalarsVarianceList():
                if s in slist: slist.remove(s)
            if slist != []:
                for s in slist:
                    diff_choice = fcm.m_sca.getScalarDiffusivityChoice(s)
                    if diff_choice == 'user_law':
                        dname = fcm.m_sca.getScalarDiffusivityName(s)
                        exp, req, sca, sym, = \
                        fcm.getFormulaComponents('scalar_diffusivity',
                                                 scalar=s)
                        self.init_block('vol', 'all_cells', dname,
                                        exp, req, sym, sca)
                        for zone in vol_zones:
                            zname = zone.getLabel()
                            z_id = str(zone.getCodeNumber())
                            if zname != "all_cells" and \
                            zone.isNatureActivated('physical_properties'):
                                exp, req, sca, sym = \
                                fcm.getFormulaComponents('scalar_diffusivity', scalar=s, zone=z_id)
                                self.init_block('vol', zname, dname,
                                                exp, req, sym, sca)

            # GroundWater Flows Law
            glm = None

            for zone in vol_zones:
                z_id = str(zone.getCodeNumber())
                zone_name = zone.getLabel()
                nature_list = zone.getNatureList()

                if "groundwater_law" in nature_list:
                    if not glm:
                        glm = GroundwaterLawModel(self.case)
                    if zone.getNature()['groundwater_law'] == 'on':
                        if glm.getGroundwaterLawModel(z_id) == 'user':
                            exp, req, sym = glm.getGroundwaterLawFormulaComponents(z_id)
                            self.init_block('vol', zone_name,
                                            'capacity+saturation+permeability',
                                            exp, req, sym, [])

            from code_saturne.model.GroundwaterModel import GroundwaterModel
            # Ground water model enabled ?
            gwm = not (GroundwaterModel(self.case).getGroundwaterModel() == 'off')

        elif self.module_name == 'neptune_cfd':
            from code_saturne.model.ThermodynamicsModel import ThermodynamicsModel
            from code_saturne.model.MainFieldsModel import MainFieldsModel
            from code_saturne.model.InterfacialEnthalpyModel import InterfacialEnthalpyModel

            tm = ThermodynamicsModel(self.case)
            mfm = MainFieldsModel(self.case)
            iem = InterfacialEnthalpyModel(self.case)

            authorized_fields = ['density', 'molecular_viscosity',
                                 'specific_heat', 'thermal_conductivity']

            compressible_fields = ['d_rho_d_P', 'd_rho_d_h']

            gas_liq_fields = ['SaturationTemperature',
                              'SaturationEnthalpyLiquid', 'SaturationEnthalpyGas',
                              'LatentHeat', 'd_Tsat_d_P',
                              'd_Hsat_d_P_Liquid', 'd_Hsat_d_P_Gas']

            user_gas_liq_fields = False
            # surface tension
            if tm:
                ## Deactivated for 7.0, might be reactivated in the future
                # if tm.getPropertyMode('none', 'surface_tension') == 'user_law':
                #     name = 'SurfaceTension'
                #     exp, req, sca, sym = tm.getFormulaComponents('none',
                #                                                  'surface_tension')
                #     self.init_block('vol', 'all_cells', name,
                #                     exp, req, sym, sca)

                for fieldId in mfm.getFieldIdList():
                    field = mfm.getFieldFromId(fieldId)
                    if tm.getMaterials(fieldId) == 'user_material':
                        for fk in authorized_fields:
                            if tm.getPropertyMode(fieldId, fk) == 'user_law':
                                name = fk + '_' + str(fieldId)
                                for zone in vol_zones:
                                    zname = zone.getLabel()
                                    z_id = str(zone.getCodeNumber())
                                    if zone.isNatureActivated('physical_properties'):
                                        exp, req, sca, sym = tm.getFormulaComponents(fieldId,fk,zone=z_id)
                                        self.init_block('vol', zname, name,
                                                        exp, req, sym, sca)

                        if field.compressible == 'on':
                            for fk in compressible_fields:
                                name = fk + '_' + str(fieldId)
                                for zone in vol_zones:
                                    zname = zone.getLabel()
                                    z_id = str(zone.getCodeNumber())
                                    if zone.isNatureActivated('physical_properties'):
                                        exp, req, sca, sym = tm.getFormulaComponents(fieldId,fk,zone=z_id)
                                        self.init_block('vol', zname, name,
                                                        exp, req, sym, sca)

                        # Temperature as a function of enthalpy
                        if field.enthalpy_model != 'off':
                            name = 'temperature_' + str(fieldId)
                            for zone in vol_zones:
                                zname = zone.getLabel()
                                z_id = str(zone.getCodeNumber())
                                if zone.isNatureActivated('physical_properties'):
                                    exp, req, sca, sym = \
                                        tm.getFormulaComponents(fieldId,
                                                                'temperature',
                                                                zone=z_id)
                                    self.init_block('vol', zname, name,
                                                    exp, req, sym, sca)

            # User properties for Water/Steam kind flows
            if tm:
                for cpl_field_ids in iem.getLiquidVaporCouples():
                    id_a = cpl_field_ids[0]
                    id_b = cpl_field_ids[1]
                    # If interfacial enthalpy model, check for user materials
                    if iem.getEnthalpyCoupleFieldId(id_a, id_b) != None:
                        if tm.getMethod(id_a) == "user_properties" and \
                                tm.getMethod(id_b) == "user_properties":
                            user_gas_liq_fields = True

            if user_gas_liq_fields:
                for fk in gas_liq_fields:
                    for zone in vol_zones:
                        zname = zone.getLabel()
                        z_id = str(zone.getCodeNumber())
                        if zone.isNatureActivated('physical_properties'):
                            exp, req, sca, sym = tm.getFormulaComponents('none', fk, zone=z_id)
                            self.init_block('vol', zname, fk,
                                            exp, req, sym, sca)

        # Porosity for both solvers
        from code_saturne.model.PorosityModel import PorosityModel

        # ALE mesh viscosity
        from code_saturne.model.MobileMeshModel import MobileMeshModel
        ale_model = MobileMeshModel(self.case)
        if ale_model.getMethod() != 'off':
            exp, req, sca, sym = ale_model.getFormulaViscComponents()
            self.init_block('vol', 'all_cells', 'mesh_viscosity',
                            exp, req, sym, sca)

        if not gwm:
            prm = PorosityModel(self.case)
            for zone in vol_zones:
                z_id = zone.getCodeNumber()
                zone_name = zone.getLabel()
                nature_list = zone.getNatureList()
                if 'porosity' in nature_list:
                    if zone.getNature()['porosity'] == 'on':
                        fname = 'porosity'
                        if prm.getPorosityModel(z_id) == 'anisotropic':
                            fname += '+tensorial_porosity'
                        exp, req, known_fields, sym = \
                        prm.getPorosityFormulaComponents(z_id)

                        self.init_block('vol', zone_name, fname,
                                        exp, req, sym, known_fields)

    # ---------------------------------------------------------------------------

    def generate_boundary_code(self):

        from code_saturne.model.NotebookModel import NotebookModel

        if self.module_name == 'code_saturne':
            from code_saturne.model.LocalizationModel import LocalizationModel
            from code_saturne.model.Boundary import Boundary
            from code_saturne.model.TurbulenceModel import TurbulenceModel

            blm = LocalizationModel('BoundaryZone', self.case)
            tm = TurbulenceModel(self.case)

            for zone in blm.getZones():

                boundary = Boundary(zone._nature, zone._label, self.case)

                # ALE: imposed mesh velocity
                c = boundary.getALEChoice()
                if c == "fixed_velocity":
                    sym = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                    for (name, val) in NotebookModel(self.case).getNotebookList():
                        sym.append((name, 'value (notebook) = ' + str(val)))
                    req = ['mesh_velocity[0]', 'mesh_velocity[1]', 'mesh_velocity[2]']
                    exp = boundary.getALEFormula()

                    name = 'mesh_velocity'
                    self.init_block('bnd', zone._label, name,
                                    exp, req, sym, known_fields=[],
                                    condition=c)
                elif c == "fixed_displacement":
                    sym = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                    for (name, val) in NotebookModel(self.case).getNotebookList():
                        sym.append((name, 'value (notebook) = ' + str(val)))
                    req = ['mesh_displacement[0]',
                           'mesh_displacement[1]',
                           'mesh_displacement[2]']
                    exp = boundary.getALEFormula()

                    name = 'mesh_velocity'
                    self.init_block('bnd', zone._label, name,
                                    exp, req, sym, known_fields=[],
                                    condition=c)

                if zone._nature == "symmetry":
                    continue

                # Velocity for inlets
                if 'inlet' in zone._nature and zone._nature != 'free_inlet_outlet':
                    c = boundary.getVelocityChoice()
                    if '_formula' in c:
                        sym = ['t', 'dt', 'iter', 'surface']
                        if c == 'norm_formula':
                            req = ['u_norm']
                            sym += ['x', 'y', 'z']
                        elif c == 'flow1_formula':
                            req = ['q_m']
                        elif c == 'flow2_formula':
                            req = ['q_v']

                        for (name, val) in NotebookModel(self.case).getNotebookList():
                            sym.append((name, 'value (notebook) = ' + str(val)))

                        name = 'velocity'

                        exp = boundary.getVelocity()
                        self.init_block('bnd', zone._label, name,
                                        exp, req, sym, known_fields=[],
                                        condition=c)

                    d = boundary.getDirectionChoice()
                    if d == 'formula':
                        req  = ['dir_x', 'dir_y', 'dir_z']
                        exp  = boundary.getDirection('direction_formula')
                        sym = ['x', 'y', 'z', 't', 'dt', 'iter']

                        for (name, val) in NotebookModel(self.case).getNotebookList():
                            sym.append((name, 'value (notebook) = ' + str(val)))

                        name = 'direction'

                        self.init_block('bnd', zone._label, name,
                                        exp, req, sym,
                                        [], condition=d)

                    # Turbulence
                    tc = boundary.getTurbulenceChoice()
                    if tc == 'formula':
                        turb_model = tm.getTurbulenceModel()
                        sym = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']

                        # Make sure the formula is "active" (i.e. is
                        # required by the current turbulence model).
                        name = None

                        for (k, val) in NotebookModel(self.case).getNotebookList():
                            sym.append((k, 'value (notebook) = ' + str(val)))

                        if turb_model in ('k-epsilon', 'k-epsilon-PL'):
                            name = 'turbulence_ke'
                            req  = ['k', 'epsilon']
                        elif turb_model in ('Rij-epsilon', 'Rij-SSG'):
                            name = 'turbulence_rije'
                            # Careful! The order of rij components must be the
                            # same as in the code (r23 before r13)
                            req  = ['r11', 'r22', 'r33',
                                    'r12', 'r23', 'r13',
                                    'epsilon']
                        elif turb_model == 'Rij-EBRSM':
                            name = 'turbulence_rij_ebrsm'
                            # Careful! The order of rij components must be the
                            # same as in the code (r23 before r13)
                            req  = ['r11', 'r22', 'r33',
                                    'r12', 'r23', 'r13',
                                    'epsilon', 'alpha']
                        elif turb_model == 'v2f-BL-v2/k':
                            name = 'turbulence_v2f'
                            req  = ['k', 'epsilon', 'phi', 'alpha']
                        elif turb_model == 'k-omega-SST':
                            name = 'turbulence_kw'
                            req  = ['k', 'omega']
                        elif turb_model == 'Spalart-Allmaras':
                            name = 'turbulence_spalart'
                            req  = ['nu_tilda']

                        if name is not None:
                            exp = boundary.getTurbFormula()
                            self.init_block('bnd', zone._label, name,
                                            exp, req, sym,
                                            [], condition=tc)

                # Specific free_inlet_outlet head loss
                if zone._nature == 'free_inlet_outlet':
                    name = "head_loss"
                    req  = ['K']
                    sym  = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                    for (nb_var, val) in NotebookModel(self.case).getNotebookList():
                        sym.append((nb_var, 'value (notebook) = ' + str(val)))

                    exp  = boundary.getHeadLossesFormula()
                    self.init_block('bnd', zone._label, name,
                                    exp, req, sym,
                                    [], condition='formula')

                # Hydraulic head for groundwater flow
                if zone._nature == 'groundwater':
                    c = boundary.getHydraulicHeadChoice()
                    sym  = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                    for (name, val) in NotebookModel(self.case).getNotebookList():
                        sym.append((name, 'value (notebook) = ' + str(val)))

                    if c == 'dirichlet_formula':
                        name = 'hydraulic_head'
                        req  = ['H']
                        exp  = boundary.getHydraulicHeadFormula()
                        self.init_block('bnd', zone._label, name,
                                        exp, req, sym,
                                        [], condition=c)

                # Scalars
                scalar_list = boundary.getDefinedScalarFormulaList()

                for e in scalar_list:
                    sca = e[0]
                    c = e[1]
                    exp = e[2]
                    sym  = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                    for (name, val) in NotebookModel(self.case).getNotebookList():
                        sym.append((name, 'value (notebook) = ' + str(val)))

                    if c == 'dirichlet_formula':
                        if sca in ('vec_potential',):
                            req = [sca+'[0]', sca+'[1]', sca+'[2]']
                        else:
                            req = [sca]
                    elif c == 'neumann_formula':
                        req = ['flux']
                    elif c == 'exchange_coefficient_formula':
                        req = [sca, 'hc']
                    else:
                        continue

                    self.init_block('bnd', zone._label, sca,
                                    exp, req, sym, [],
                                    condition=c)

        else:
            from code_saturne.model.LocalizationModel import LocalizationModel
            from code_saturne.model.BoundaryNeptune import Boundary
            from code_saturne.model.MainFieldsModel import MainFieldsModel
            from code_saturne.model.TurbulenceNeptuneModel import TurbulenceModel
            from code_saturne.model.SpeciesModel import SpeciesModel

            blm = LocalizationModel("BoundaryZone", self.case)
            mfm = MainFieldsModel(self.case)
            tm  = TurbulenceModel(self.case)
            spm = SpeciesModel(self.case)

            for zone in blm.getZones():
                if "inlet" in zone.getNature():
                    for fId in mfm.getFieldIdList():
                        boundary = Boundary(zone.getNature(),
                                            zone.getLabel(),
                                            self.case,
                                            fId)

                        # Velocity
                        c = boundary.getVelocityChoice(fId)
                        if '_formula' in c:
                            sym = []
                            if c == 'norm_formula':
                                req = ['u_norm']
                                sym.extend(('x', 'y', 'z'))
                            elif c == 'flow1_formula':
                                req = ['q_m']
                            sym.extend(('t', 'dt', 'iter', 'surface'))

                            for (name, val) in NotebookModel(self.case).getNotebookList():
                                sym.append((name, 'value (notebook) = ' + str(val)))

                            exp = boundary.getVelocity(fId)

                            self.init_block('bnd',
                                            zone.getLabel(),
                                            'velocity_'+str(fId),
                                            exp,
                                            req,
                                            sym,
                                            [],
                                            condition=c)

                        # Velocity direction
                        d = boundary.getDirectionChoice(fId)
                        if d == 'formula':
                            exp = boundary.getDirection(fId, 'direction_formula')
                            req = ['dir_x', 'dir_y', 'dir_z']
                            sym = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                            for (name, val) in NotebookModel(self.case).getNotebookList():
                                sym.append((name, 'value (notebook) = ' + str(val)))

                            self.init_block('bnd',
                                            zone.getLabel(),
                                            'direction_'+str(fId),
                                            exp,
                                            req,
                                            sym,
                                            [],
                                            condition = d)

                        # Turbulence
                        tc = boundary.getTurbulenceChoice(fId)
                        turb_model = tm.getTurbulenceModel(fId)
                        if tc == 'formula' and turb_model not in ['none', 'q2-q12-tchen', 'r2-r12-tchen']:
                            exp, reqo, sym = boundary.getTurbFormulaComponents(fId,
                                                                              turb_model)
                            if turb_model in('k-epsilon', 'k-epsilon_linear_production'):
                                name = 'turbulence_ke_%s' % (fId)
                            elif turb_model in ('k-omega-SST'):
                                name = 'turbulence_kw_%s' % (fId)
                            elif turb_model in ('rij-epsilon_ssg', 'rij-epsilon_ebrsm'):
                                name = 'turbulence_rije_%s' % (fId)
                            elif turb_model in ('tchen', 'q2-q12'):
                                name = 'turbulence_tchen_%s' % (fId)
                            elif turb_model in ('r2-q12'):
                                name = 'turbulence_r2q12_%s' % (fId)
                            elif turb_model in ('r2-r12-tchen'):
                                name = 'turbulence_r2r12_%s' % (fId)

                            if type(reqo[0]) == tuple:
                                req = [r[0] for r in reqo]
                            else:
                                req = reqo

                            self.init_block('bnd',
                                            zone.getLabel(),
                                            name,
                                            exp,
                                            req,
                                            sym,
                                            [],
                                            condition=tc)
                # Thermal conditions
                if zone.getNature() in ['inlet', 'outlet']:
                    for fId in mfm.getFieldIdList():
                        boundary = Boundary(zone.getNature(),
                                            zone.getLabel(),
                                            self.case,
                                            fId)

                        c = boundary.getEnthalpyChoice(fId)
                        if '_formula' in c:
                            sym  = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                            if c == "dirichlet_formula":
                                req = ['enthalpy']
                            elif c == "flux_formula":
                                req = ['flux']
                            elif c in ['temperature_formula', 'timp_K_formula']:
                                req = ['temperature']

                            for (name, val) in NotebookModel(self.case).getNotebookList():
                                sym.append((name, 'value (notebook) = ' + str(val)))

                            exp = boundary.getEnthalpy(fId)
                            self.init_block('bnd',
                                            zone.getLabel(),
                                            'enthalpy_'+str(fId),
                                            exp,
                                            req,
                                            sym,
                                            [],
                                            condition=c)

                if zone.getNature() == "wall":
                    boundary = Boundary(zone.getNature(),
                                        zone.getLabel(),
                                        self.case,
                                        'none')

                    # ALE: imposed mesh velocity
                    c = boundary.getALEChoice()

                    if c == "fixed_velocity":
                        sym = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                        for (name, val) in NotebookModel(self.case).getNotebookList():
                            sym.append((name, 'value (notebook) = ' + str(val)))
                        req = ['mesh_velocity[0]', 'mesh_velocity[1]', 'mesh_velocity[2]']
                        exp = boundary.getALEFormula()

                        name = 'mesh_velocity'
                        self.init_block('bnd', zone._label, name,
                                        exp, req, sym, known_fields=[],
                                        condition=c)
                    elif c == "fixed_displacement":
                        sym = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                        for (name, val) in NotebookModel(self.case).getNotebookList():
                            sym.append((name, 'value (notebook) = ' + str(val)))
                        req = ['mesh_displacement[0]',
                               'mesh_displacement[1]',
                               'mesh_displacement[2]']
                        exp = boundary.getALEFormula()

                        name = 'mesh_velocity'
                        self.init_block('bnd', zone._label, name,
                                        exp, req, sym, known_fields=[],
                                        condition=c)

                    c = boundary.getEnthalpyChoice('none')
                    if '_formula' in c:
                        sym  = ['x', 'y', 'z', 't', 'dt', 'iter', 'surface']
                        if c == "dirichlet_formula":
                            req = ['enthalpy']
                        elif c == "flux_formula":
                            req = ['flux']
                        elif c in ['temperature_formula', 'timp_K_formula']:
                            req = ['temperature']

                        for (name, val) in NotebookModel(self.case).getNotebookList():
                            sym.append((name, 'value (notebook) = ' + str(val)))

                        exp = boundary.getEnthalpy('none')
                        self.init_block('bnd',
                                        zone.getLabel(),
                                        'enthalpy',
                                        exp,
                                        req,
                                        sym,
                                        [],
                                        condition=c)

                # Scalars
                for sp_id in mfm.getFieldIdList(include_none=True):
                    boundary = Boundary(zone.getNature(),
                                        zone.getLabel(),
                                        self.case,
                                        sp_id)
                    if boundary.getNature() != "symmetry":
                        for _s in spm.getScalarByFieldId(sp_id):
                            c = boundary.getScalarChoice(sp_id, _s)
                            if 'formula' == c[-7:]:
                                exp, req, sym = \
                                    boundary.getScalarFormulaComponents(sp_id, _s)
                                self.init_block('bnd',
                                                zone.getLabel(),
                                                spm.getScalarLabelByName(_s),
                                                exp,
                                                req,
                                                sym,
                                                [],
                                                condition=c)

    # ---------------------------------------------------------------------------

    def generate_source_terms_code(self, vol_zones=None):

        if vol_zones == None:
            from code_saturne.model.LocalizationModel import LocalizationModel
            vlm = LocalizationModel('VolumicZone', self.case)
            vol_zones = vlm.getZones()

        if self.module_name == 'code_saturne':
            from code_saturne.model.SourceTermsModel import SourceTermsModel
            from code_saturne.model.GroundwaterModel import GroundwaterModel
            from code_saturne.model.DefineUserScalarsModel import DefineUserScalarsModel

            stm = SourceTermsModel(self.case)
            gwm = GroundwaterModel(self.case)

            for zone in vol_zones:
                z_id = str(zone.getCodeNumber())
                zone_name = zone.getLabel()

                nature_list = zone.getNatureList()

                if 'momentum_source_term' in nature_list:
                    if zone.getNature()['momentum_source_term'] == 'on':
                        if gwm.getGroundwaterModel() == 'off':
                            exp, req, sym = stm.getMomentumFormulaComponents(z_id)
                            knf = [('rho','density')]
                            self.init_block('src', zone_name, "momentum",
                                            exp, req, sym, knf,
                                            source_type="momentum_source_term")
                        else:
                            exp, req, sym = stm.getRichardsFormulaComponents(z_id)
                            self.init_block('src', zone_name, 'richards',
                                            exp, req, sym, [],
                                            source_type="momentum_source_term")

                if 'scalar_source_term' in nature_list:
                    if zone.getNature()['scalar_source_term'] == 'on':
                        sca_list = DefineUserScalarsModel(self.case).getUserScalarNameList()
                        if gwm.getGroundwaterModel() == 'off':
                            for sca in sca_list:
                                exp, req, sym, knf = stm.getSpeciesFormulaComponents(z_id, sca)

                                self.init_block('src', zone_name, sca,
                                                exp, req, sym, knf,
                                                source_type="scalar_source_term")
                        else:
                            for sca in sca_list:
                                exp, req, sym, knf = \
                                stm.getGroundWaterSpeciesFormulaComponents(z_id, sca)

                                self.init_block('src', zone_name, sca,
                                                exp, req, sym, knf,
                                                source_type="scalar_source_term")

                if 'thermal_source_term' in nature_list:
                    if zone.getNature()['thermal_source_term'] == 'on':
                        th_sca_name = stm.therm.getThermalScalarName()
                        exp, req, sym, knf = stm.getThermalFormulaComponents(z_id,
                                                                             th_sca_name)

                        self.init_block('src', zone_name, th_sca_name,
                                        exp, req, sym, knf,
                                        source_type="thermal_source_term")
        else:
            from code_saturne.model.LocalizationModel import LocalizationModel
            from code_saturne.model.MainFieldsSourceTermsModel import MainFieldsSourceTermsModel

            stm = MainFieldsSourceTermsModel(self.case)

            for zone in vol_zones:
                z_id = str(zone.getCodeNumber())
                zone_name = zone.getLabel()

                nature_list = zone.getNatureList()
                if 'thermal_source_term' in nature_list:
                    if zone.getNature()['thermal_source_term'] == 'on':
                        for fId in stm.mfm.getFieldIdList():
                            exp, req, sym, knf = stm.getThermalFormulaComponents(z_id,
                                                                                 fId,
                                                                                 'enthalpy')
                            self.init_block('src', zone_name,
                                            'enthalpy_'+str(fId),
                                            exp, req, sym, knf,
                                            source_type='thermal_source_term')

    # ---------------------------------------------------------------------------

    def generate_initialize_code(self, vol_zones=None):

        if vol_zones == None:
            from code_saturne.model.LocalizationModel import LocalizationModel
            vlm = LocalizationModel('VolumicZone', self.case)
            vol_zones = vlm.getZones()

        if self.module_name == 'code_saturne':
            from code_saturne.model.InitializationModel import InitializationModel
            from code_saturne.model.CompressibleModel import CompressibleModel
            from code_saturne.model.DefineUserScalarsModel import DefineUserScalarsModel
            im = InitializationModel(self.case)
            cpm = CompressibleModel(self.case)

            for zone in vol_zones:
                if zone.getNature()['initialization'] == 'on':
                    z_id = str(zone.getCodeNumber())
                    zone_name = zone.getLabel()

                    # Velocity
                    if im.getVelocityFormula(z_id):
                        exp, req, sym = im.getVelocityFormulaComponents(z_id)
                        self.init_block('ini', zone_name, 'velocity',
                                        exp, req, sym, [])

                    # Void fraction
                    if im.getVoidFractionFormula(z_id):
                        exp, req, sym = im.getVoidFractionFormulaComponents(z_id)
                        self.init_block('ini', zone_name, 'void_fraction',
                                        exp, req, sym, [])

                    # Turbulence
                    tin = im.node_turb.xmlGetNode('initialization', zone_id=z_id)
                    if tin:
                        if tin['choice'] == 'formula':
                            tmodel = im.node_turb['model']
                            exp, req, sym = im.getTurbFormulaComponents(z_id, tmodel)
                            self.init_block('ini', zone_name, 'turbulence',
                                            exp, req, sym, [])

                    # Thermal
                    node_t = im.node_scalartherm.xmlGetNode('variable')
                    if node_t:
                        th_formula = im.getThermalFormula(z_id)
                        if th_formula:
                            exp, req, sym, knf = im.getThermalFormulaComponents(z_id)
                            self.init_block('ini', zone_name, 'thermal',
                                            exp, req, sym, knf)

                    # HydraulicHead
                    if im.node_veloce.xmlGetNode('variable', name = 'hydraulic_head'):
                        if im.getHydraulicHeadFormula(z_id):
                            exp, req, sym = im.getHydraulicHeadFormulaComponents(z_id)
                            self.init_block('ini', zone_name, 'hydraulic_head',
                                            exp, req, sym, [])

                    if cpm.getCompressibleModel() != 'off':
                        # Pressure
                        if im.getPressureStatus(z_id) != 'off':
                            exp, req, sym = im.getPressureFormulaComponents(z_id)
                            self.init_block('ini', zone_name, 'pressure',
                                            exp, req, sym, [])

                        # Density
                        if im.getDensityStatus(z_id) != 'off':
                            exp, req, sym = im.getDensityFormulaComponents(z_id)
                            self.init_block('ini', zone_name, 'density',
                                            exp, req, sym, [])

                        # Temperature
                        if im.getTemperatureStatus(z_id) != 'off':
                            exp, req, sym = im.getTemperatureFormulaComponents(z_id)
                            self.init_block('ini', zone_name, 'temperature',
                                            exp, req ,sym, [])

                        # Energy
                        if im.getEnergyStatus(z_id) != 'off':
                            exp, req, sym = im.getEnergyFormulaComponents(z_id)
                            self.init_block('ini', zone_name, 'energy',
                                            exp, req, sym, [])

                    # Species
                    usm = DefineUserScalarsModel(self.case)
                    for scalar in usm.getUserScalarNameList():
                        if im.getSpeciesFormula(z_id, scalar):
                            exp, req, sym, knf = im.getSpeciesFormulaComponents(z_id, scalar)
                            self.init_block('ini', zone_name, scalar,
                                            exp, req, sym, knf)

                    # Meteo
                    node_atmo = im.models.xmlGetNode('atmospheric_flows')
                    if node_atmo:
                        for mscalar in usm.getMeteoScalarsNameList():
                            node = node_atmo.xmlGetNode('variable', name=mscalar)
                            if node.xmlGetString('formula', zone_id=z_id):
                                exp,req,sym = im.getMeteoFormulaComponents(z_id,mscalar)
                                self.init_block('ini', zone_name, mscalar,
                                                exp, req, sym, [])

                    # Combustion
                    node_gas = im.models.xmlGetNode('gas_combustion')
                    if node_gas:
                        for mscalar in usm.getGasCombScalarsNameList():
                            node = node_gas.xmlGetNode('variable', name=mscalar)
                            if node.xmlGetString('formula', zone_id=z_id):
                                exp,req,sym = im.getCombustionFormulaComponents(z_id,mscalar)
                                self.init_block('ini', zone_name, mscalar,
                                                exp, req, sym, [])

        else:
            from code_saturne.model.MainFieldsModel import MainFieldsModel
            from code_saturne.model.MainFieldsInitializationModel import MainFieldsInitializationModel
            from code_saturne.model.NonCondensableModel import NonCondensableModel
            from code_saturne.model.SpeciesModel import SpeciesModel
            from code_saturne.model.TurbulenceNeptuneModel import TurbulenceModelsDescription
            from code_saturne.model.TurbulenceNeptuneModel import TurbulenceModel

            mfm = MainFieldsModel(self.case)
            mfi = MainFieldsInitializationModel(self.case)
            ncm = NonCondensableModel(self.case)
            spm = SpeciesModel(self.case)

            for zone in vol_zones:
                if zone.getNature()['initialization'] == 'on':
                    z_id = str(zone.getCodeNumber())
                    zone_name = zone.getLabel()

                    # Pressure
                    exp, req, sym = mfi.getPressureFormulaComponents(z_id)
                    self.init_block('ini', zone_name,
                                    'pressure',
                                    exp, req, sym, [])

                    for fId in mfm.getFieldIdList():

                        # Velocity
                        exp, req, sym = mfi.getFormulaComponents(z_id,
                                                                 fId,
                                                                 'velocity')
                        self.init_block('ini', zone_name, 'velocity_'+str(fId),
                                        exp, req, sym, [])

                        # Volume fraction
                        exp, req, sym = mfi.getFormulaComponents(z_id,
                                                                 fId,
                                                                 'volume_fraction')
                        self.init_block('ini', zone_name,
                                        'volume_fraction_'+str(fId),
                                        exp, req, sym, [])

                        # Turbulence get turbModel here before
                        turbModel = \
                                TurbulenceModel(self.case).getTurbulenceModel(fId)
                        if turbModel != 'none' and 'mixing' not in turbModel:
                            exp, req, sym = \
                                    TurbulenceModel(self.case).getFormulaComponents(z_id,
                                                                                    fId,
                                                                                    turbModel)

                            # Play with variablesDim to find the start and end index for the
                            # required variables
                            self.init_block('ini', zone_name,
                                            'turbulence_'+str(fId),
                                            exp,
                                            req,
                                            sym,
                                            [])

                        # Enthalpy (only if energy resolution is activated)
                        if mfm.getFieldFromId(fId).enthalpy_model != 'off':
                            if mfi.getEnergyModel(z_id, fId) != 'hsat_P':
                                exp, req, sym = mfi.getFormulaComponents(z_id,
                                                                         fId,
                                                                         'enthalpy')
                                self.init_block('ini', zone_name,
                                                'enthalpy_'+str(fId),
                                                exp, req, sym, [])

                        # Non condensables
                        for nc in ncm.getNonCondensableByFieldId(fId):
                            exp, req, sym = \
                                mfi.getNonCondensableFormulaComponents(z_id,
                                                                       fId,
                                                                       nc)
                            self.init_block('ini', zone_name,
                                            ncm.getNonCondLabel(nc),
                                            exp, req, sym, [])

                    # Species are treated apart since they can be phase-independent
                    for fId in mfm.getFieldIdList(include_none=True):
                        for s in spm.getScalarByFieldId(fId):
                            exp, req, sym = \
                                mfi.getScalarFormulaComponents(z_id,
                                                               fId,
                                                               s)
                            self.init_block('ini', zone_name,
                                            spm.getScalarLabelByName(s),
                                            exp, req, sym, [])

    # ---------------------------------------------------------------------------

    def generate_immersed_boundaries_code(self):

        from code_saturne.model.ImmersedBoundariesModel import ImmersedBoundariesModel
        ibm = ImmersedBoundariesModel(self.case)

        if ibm.getOnOff() == 'on':
            for objId in range(1,ibm.getNumberOfObjects()+1):

                object_name = ibm.getObjectName(objId)

                if ibm.getObjectMethod(objId) == 'defined by function':
                    exp, req, sym = ibm.getIBMFormulaComponents(objId-1)
                    self.init_block('ibm', object_name, 'porosity',
                                    exp, req, sym, [])

                #inertia (mass is constant)
                if ibm.getObjectModelingMode(objId) == 'mass':
                    exp, req, sym = ibm.getFormulaInertia(objId-1)
                    self.init_block('ibm_fsi', object_name, 'porous_inertia',
                                    exp, req, sym, [])

                #rho
                if (ibm.getObjectCHT(objId) == "on" or ibm.getObjectFSI(objId) == "on"):
                    if ibm.getObjectPropertyMode(objId, 'density') == 'user_law':
                        exp, req, sym = ibm.getFormulaRho(objId-1)
                        self.init_block('ibm_vol', object_name, 'porous_density',
                                    exp, req, sym, [])

                #CHT
                if (ibm.getObjectCHT(objId) == "on"):
                    #CP
                    if ibm.getObjectPropertyMode(objId, 'specific_heat') == 'user_law':
                        exp, req, sym = ibm.getFormulaCp(objId-1)
                        self.init_block('ibm_vol', object_name, 'porous_specific_heat',
                                        exp, req, sym, [])

                    #lambda
                    if ibm.getObjectPropertyMode(objId, 'thermal_conductivity') == 'user_law':
                        exp, req, sym = ibm.getFormulaAl(objId-1)
                        self.init_block('ibm_vol', object_name, 'porous_thermal_conductivity',
                                        exp, req, sym, [])

                #Thermal source term (only if CHT)
                if ibm.getObjectThermalSourceTerm(objId) == 'on':
                    exp, req, sym = ibm.getFormulaThermalSourceTerm(objId-1)
                    self.init_block('ibm_vol', object_name, 'porous_thermal_st',
                                    exp, req, sym, [])

                # Moving
                if ibm.getObjectMoving(objId) == 'imposed':
                    exp, req, sym = ibm.getImposedMovingFormulaComponents(objId-1)
                    self.init_block('ibm_vol', object_name, 'porous_velocity',
                                    exp, req, sym, [])

                # initial temperature
                if (ibm.getObjectInit(objId) == 'on'):
                    exp, req, sym = ibm.getFormulaTemperature(objId-1)
                    self.init_block('ibm_vol', object_name, 'porous_temperature',
                                    exp, req, sym, [])

                # Boundary
                energy_mode = ibm.getObjectBoundaryEnergyMode(objId)
                if (energy_mode in ['temperature_formula', 'flux_formula']):
                    exp, req, sym = ibm.getFormulaBoundaryEnergy(objId-1, energy_mode)
                    if (energy_mode == 'temperature_formula'):
                        self.init_block('ibm_vol', object_name, 'boundary_temperature',
                                        exp, req, sym, [])
                    elif (energy_mode == 'flux_formula'):
                        self.init_block('ibm_vol', object_name, 'boundary_heat_flux',
                                        exp, req, sym, [])

    # ---------------------------------------------------------------------------

    def generate_fsi_ic_code(self):
        # ALE enabled ?

        from code_saturne.model.NotebookModel import NotebookModel

        from code_saturne.model.LocalizationModel import LocalizationModel
        from code_saturne.model.Boundary import Boundary
        from code_saturne.model.TurbulenceModel import TurbulenceModel

        blm = LocalizationModel('BoundaryZone', self.case)

        for zone in blm.getZones():

            boundary = Boundary(zone._nature, zone._label, self.case)

            # ALE: imposed mesh velocity
            c = boundary.getALEChoice()
            if c != "internal_coupling":
                continue

            boundary = Boundary("coupling_mobile_boundary", zone._label, self.case)

            m_mat = boundary.getMassMatrix()
            k_mat = boundary.getStiffnessMatrix()
            c_mat = boundary.getDampingMatrix()
            f_force = boundary.getFluidForce()

            sym = ['t', 'dt', 'iter']
            for (name, val) in NotebookModel(self.case).getNotebookList():
                sym.append((name, 'value (notebook) = ' + str(val)))

            req = ['m11', 'm12', 'm13', 'm21', 'm22', 'm23', 'm31', 'm32', 'm33']
            self.init_block('fsi', zone._label, 'mass_matrix',
                            m_mat, req, sym, known_fields=[])

            req = ['k11', 'k12', 'k13', 'k21', 'k22', 'k23', 'k31', 'k32', 'k33']
            self.init_block('fsi', zone._label, 'stiffness_matrix',
                            k_mat, req, sym, known_fields=[])

            req = ['c11', 'c12', 'c13', 'c21', 'c22', 'c23', 'c31', 'c32', 'c33']
            self.init_block('fsi', zone._label, 'damping_matrix',
                            c_mat, req, sym, known_fields=[])

            sym = ['t', 'dt', 'iter', 'fluid_fx', 'fluid_fy', 'fluid_fz']
            for (name, val) in NotebookModel(self.case).getNotebookList():
                sym.append((name, 'value (notebook) = ' + str(val)))
            req = ['fx', 'fy', 'fz']

            self.init_block('fsi', zone._label, 'fluid_force',
                            f_force, req, sym, known_fields=[])

    # ---------------------------------------------------------------------------

    def generate_post_profile_code(self):
        # Output writer activation

        from code_saturne.model.NotebookModel import NotebookModel
        from code_saturne.model.ProfilesModel import ProfilesModel

        pfm = ProfilesModel(self.case)

        for l in pfm.getProfilesLabelsList():

            formula = pfm.getFormula(l)

            if not formula:
                continue

            npts = pfm.getNbPoint(l)

            sym = ['s']
            for (name, val) in NotebookModel(self.case).getNotebookList():
                sym.append((name, 'value (notebook) = ' + str(val)))

            req = ['x', 'y', 'z']

            self.init_block('pfl', l, 'coordinates',
                            formula, req, sym, known_fields=[])

    # ---------------------------------------------------------------------------

    def generate_writer_activation_code(self):
        # Output writer activation

        from code_saturne.model.NotebookModel import NotebookModel
        from code_saturne.model.OutputControlModel import OutputControlModel

        ocm = OutputControlModel(self.case)

        for writer_id in ocm.getWriterIdList():

            formula = None
            frequency_choice = ocm.getWriterFrequencyChoice(writer_id)
            if frequency_choice == "formula":
                formula = ocm.getWriterFrequency(writer_id)

            if not formula:
                continue

            sym = ['t', 'iter']
            for (name, val) in NotebookModel(self.case).getNotebookList():
                sym.append((name, 'value (notebook) = ' + str(val)))

            req = ['is_active']

            self.init_block('pwa', str(writer_id), 'activation',
                            formula, req, sym, known_fields=[])

    def generate_post_calculator_code(self):
        # Postprocessing calculator

        from code_saturne.model.UserCalculatorModel import UserCalculatorModel
        calculator = UserCalculatorModel(self.case)

        for name in calculator.getFunctionsNamesList():
            exp, req, known_fields, sym = calculator.getFunctionFormulaComponents(name)
            self.init_block('pca', 'None', name, exp, req, sym, known_fields)

    # ---------------------------------------------------------------------------

    def check_meg_code_syntax(self, function_name):

        if not os.path.exists(self.tmp_path):
            os.makedirs(self.tmp_path)

        if function_name in ('vol', 'bnd', 'src', 'ini', 'ibm', 'ibm_vol',
                             'ibm_fsi', 'fsi', 'pfl', 'pwa'):
            try:
                self.save_function(func_type=function_name,
                                   hard_path=self.tmp_path)
            except Exception as e:
                return -1, str(e), 1

        from code_saturne.base import cs_compile

        cwd = os.getcwd()
        os.chdir(self.tmp_path)

        out = open('comp.out', 'w')
        err = open('comp.err', 'w')

        solver = "cs_solver" + self.case['package'].config.exeext
        if self.case.module_name() == 'neptune_cfd':
            solver = "nc_solver" + self.case['package'].config.exeext

        compilation_test = cs_compile.compile_and_link(self.case['package'],
                                                       solver,
                                                       self.tmp_path,
                                                       opt_cflags='-w',
                                                       stdout=out,
                                                       stderr=err)
        out.close()
        err.close()

        n_errors = 0
        msg = ''
        if compilation_test != 0:
            errors = open('comp.err', 'r').readlines()
            for i in range(len(errors)):
                if ': ' in errors[i]:
                    msg += errors[i].split(': ')[-1].strip()+'\n'
                    n_errors += 1
            if n_errors == 0: # in case we cannot parse the output correctly
                n_errors += 1
                for i in range(len(errors)):
                    msg += errors[i].strip()+'\n'

        os.chdir(cwd)

        return compilation_test, msg, n_errors

    # ---------------------------------------------------------------------------

    def clean_tmp_dir(self):

        if os.path.exists(self.tmp_path):
            fl = os.listdir(self.tmp_path)
            for f in fl:
                self.delete_file(f, self.tmp_path)
            os.rmdir(self.tmp_path)

    # ---------------------------------------------------------------------------

    def has_meg_code(self):

        retcode = False

        if getRunType(self.case) == 'standard':
            for func_type in self.funcs.keys():
                if len(self.funcs[func_type].keys()) > 0:
                    retcode = True
                    break

        return retcode

    # ---------------------------------------------------------------------------

    def __file_path__(self, c_file_name, hard_path=None):

        # Path based on call options
        if hard_path != None:
            fpath = os.path.join(hard_path, c_file_name)
        elif self.wdir:
            fpath = os.path.join(self.wdir, c_file_name)
        else:
            fpath = os.path.join(self.case['case_path'], 'src', c_file_name)

        return fpath

    # ---------------------------------------------------------------------------

    def delete_file(self, c_file_name, hard_path=None):

        # Copy function file if needed
        fpath = self.__file_path__(c_file_name, hard_path=hard_path)

        if os.path.isfile(fpath):
            os.remove(fpath)

    # ---------------------------------------------------------------------------

    def clean_lines(self, code_to_write):

        lines = code_to_write.split('\n')
        code_to_write = ""

        for i, l in enumerate(lines):
            lines[i] = l.rstrip()

        return '\n'.join(lines)

    # ---------------------------------------------------------------------------

    def save_file(self, c_file_name, code_to_write, hard_path=None):

        if code_to_write != '':
            # Try and write the function in the src if in RESU folder
            # For debugging purposes
            try:
                fpath = self.__file_path__(c_file_name, hard_path=hard_path)
                new_file = open(fpath, 'w')
                new_file.write(self.clean_lines(code_to_write))
                new_file.close()
                return 1

            except:
                # Cant save the function. xml file will still be saved
                return 2

        # Return 0 if nothing is written for robustness
        else:
            return 0

    # ---------------------------------------------------------------------------

    def save_function(self, func_type, hard_path = None):

        # Temporary modification to handle specific cases (IBM):
        if func_type in ('ibm', 'ibm_vol', 'ibm_fsi'):
            return self.save_functions_file(func_type, hard_path)

        # Delete previous existing file
        file2write = _function_names[func_type]
        self.delete_file(file2write)

        # Check if it is a standard computation
        if getRunType(self.case) != 'standard':
            return 0

        # Generate the functions code if needed
        code_to_write = ''
        if len(self.funcs[func_type].keys()) > 0:
            code_to_write = _file_header
            #            if self.module_name != "code_saturne":
            #                code_to_write += _file_header2
            code_to_write += _file_header3
            code_to_write += _function_header[func_type]
            k_count = 0
            for key in self.funcs[func_type].keys():
                w_block = self.write_block(func_type, key)
                if w_block is None:
                    continue
                zone_name, var_name = key.split('::')
                var_name = var_name.replace("+", ", ")
                if zone_name != 'None':
                    m1 = _block_comments[func_type] % (var_name, zone_name)
                else:
                    m1 = _block_comments[func_type] % (var_name)
                m2 = '  -' + '-'*len(m1) + ' */\n\n'
                m1 = '/* ' + m1 + '\n'

                if k_count > 0:
                    code_to_write += '\n'
                code_to_write += '  ' + m1
                code_to_write += '  ' + m2
                code_to_write += w_block

                k_count += 1

            code_to_write += _file_footer

        # Write the C file if necessary
        save_status = self.save_file(file2write,
                                     code_to_write,
                                     hard_path = hard_path)

        return save_status

    # ---------------------------------------------------------------------------

    def save_functions_file(self, func_type, hard_path = None):

        # Delete previous existing
        file2write = _function_names[func_type]
        self.delete_file(file2write)

        # Check if it is a standard computation
        if getRunType(self.case) != 'standard':
            return 0

        # Generate the functions code if needed
        code_to_write = ''

        if len(self.funcs[func_type].keys()) > 0:
            code_to_write = _file_header
            code_to_write += _file_header3

            # We build the file while first writing the static internal
            # functions, then the public function (getter)
            code_to_write += _private_func_header

            getter_func_code = _function_getter_header[func_type]
            getter_func_code += '  if (object_name == nullptr)\n'
            getter_func_code += '    return nullptr;\n'
            k_count = 0

            for key in self.funcs[func_type].keys():
                w_block = self.write_block(func_type, key)
                if w_block == None:
                    continue

                objname = key.split('::')[0]
                var_name = key.split('::')[1]

                suffix = objname
                if (func_type in ["ibm_vol", "ibm_fsi"]):
                    suffix = '_'.join([objname, var_name])

                _fname = build_internal_function_header(func_type, suffix)
                code_to_write += \
                        _internal_function_prefix[func_type].format(_name=_fname,_sep_open='{')
                code_to_write += w_block
                code_to_write += _function_close

                if (func_type in ["ibm"]):
                    getter_func_code += '  else if (strcmp(object_name, "{}") == 0)\n'.format(objname)
                elif (func_type in ["ibm_vol", "ibm_fsi"]):
                    getter_func_code += (
                        '  else if (   strcmp(object_name, "{}") == 0\n'
                        '           && strcmp(var_name, "{}") == 0)\n'
                    ).format(objname, var_name)

                getter_func_code += '    return {};\n'.format(_fname)

                k_count += 1

            code_to_write += _public_func_header

            code_to_write += getter_func_code
#            code_to_write += _function_close

            code_to_write += _file_footer

        # Write the C file if necessary
        save_status = self.save_file(file2write,
                                     code_to_write,
                                     hard_path = hard_path)

        return save_status


    #---------------------------------------------------------------------------

    def save_all_functions(self):

        save_status = 0

        is_empty    = 0
        empty_exps  = []
        for func_type in self.funcs.keys():
            try:
                state = self.save_function(func_type)
            except Exception as e:
                msg = "Error while generating \"%s\" formulae.\n" % \
                        _func_short_to_long[func_type]
                msg+=str(e)
                _error_and_exit(msg)

            if state != 0:
                save_status = state

            for ek in self.funcs[func_type].keys():
                if self.funcs[func_type][ek]['exp'] in [None, ""]:
                    is_empty = 1

                    empty_exps.append({})
                    empty_exps[-1]['zone'] = ek.split('::')[0]
                    empty_exps[-1]['var']  = ek.split('::')[1]
                    empty_exps[-1]['func'] = _func_short_to_long[func_type] + ' formula'

        if is_empty == 1:
            save_status = -1

        ret = {'state':save_status,
               'exps':empty_exps,
               'nexps':len(empty_exps)}

        return ret

# -------------------------------------------------------------------------------
# End
# -------------------------------------------------------------------------------
